Observations from the station of Theeuwes

wd<-"D:/uhimax/"

uhi<-readRDS(paste0(wd,"Rdata/uhimax_obs_calc.rds"))
uhi<-uhi[which(uhi$stn %in% "UTRECHT23"),]

summary(uhi$meteo)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.4471  1.8583  2.5031  2.4961  3.1021  5.0066

Gridded data from SVF and vegetation fraction

Plotting all the data without excluding values. Extent of the area is cropped to the city center of Utrecht.

city_center<-extent(c(5.109,5.13,52.08,52.098))
city_center<-raster(city_center)
values(city_center)<-1
crs(city_center)<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
# city_center<-projectRaster(city_center,crs=crs(grd.fveg))

grd.fveg<-stack(paste0(wd,"Grids_veg_svf/greenness_utrecht_smooth_500m.grd"))

svf_utrecht<-list.files(paste0(wd,"SVF_new_subset"),pattern=".grd",full.names = TRUE)
raster_utrecht<-lapply(svf_utrecht,raster)

# mapview(raster_utrecht[[1]])

city_center<- projectRaster(city_center,crs=crs(raster_utrecht[[1]]))
grd.fveg<- projectRaster(grd.fveg,crs=crs(raster_utrecht[[1]]))

raster_utrecht_center<-mapply(crop,raster_utrecht,MoreArgs=list(y=extent(city_center)))
grd.fveg_center<-crop(grd.fveg,extent(raster_utrecht_center[[1]]))

mapview(raster_utrecht_center[[1]]) + grd.fveg_center + raster_utrecht_center[[2]] + raster_utrecht_center[[3]]
## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 3631200 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  3631200 '

Calculating the UHI according to Theeuwes

Low and high svf and vegetation fration values are excluded.

uhi_center_utrecht<-function(skyview){
  skyview<-raster::resample(skyview,grd.fveg_center,method="bilinear")
  grd.fveg_center<-crop(grd.fveg_center,extent(skyview))
  skyview<-crop(skyview,extent(grd.fveg_center))
  
  values(skyview)[values(skyview)<0.2 | values(skyview)>0.9] = NA
  values(grd.fveg_center)[values(grd.fveg_center)<0 | values(grd.fveg_center)>0.4] = NA

  uhi_center<-(2-grd.fveg_center-skyview)*median(uhi$meteo)
  return(list("svf"=skyview,"fveg"=grd.fveg_center,"meteo"=median(uhi$meteo),"uhi"=uhi_center))
  }
uhi_center<-lapply(raster_utrecht_center,uhi_center_utrecht)
mapview(uhi_center[[1]]$uhi) + uhi_center[[2]]$uhi + uhi_center[[3]]$uhi
mapview(uhi_center[[1]]$fveg)